Zmiany złożoności przestrzennej obszarów leśnych w latach 1992-2020
Wprowadzenie
Poniższy raport obejmuje przygotowanie danych oraz przebieg analizy i wyniki przygotowane dla jednej z trzech klas pokrycia terenu - obszarów leśnych. Posłużono się jedną, przykładową klasą, ze względu na skupieniu się na przedstawieniu metodyki oraz prezentacji głównych wyników. Celem analizy jest określenie zmian złożoności przestrzennej obszarów leśnych na świecie w latach 1992-2020.
Strukturę przestrzenną lokalnego krajobrazu zaprezentowanego w formie kategorycznego rastra można określić m. in. za pomocą macierzy współwystępowania, czyli relacji pomiędzy sąsiadującymi komórkami rastra. Na ich podstawie można obliczyć dwie główne metryki: entropię brzegową oraz względną informację wzajemną. Określają one kolejno różnorodność (złożoność) oraz konfigurację (jednolitość, zbitość) lokalnego krajobrazu. Czym większa wartość entropii brzegowej, tym większa różnorodność krajobrazu, czym większa wartość względnej informacji wzajemnej - tym krajobraz jest bardziej jednolity przestrzennie.
Przygotowanie danych
1.1 Wczytanie danych rastrowych
Rastry obejmujące globalne pokrycie terenu na lata 1992-2020 pozyskano z projektu ESA Climate Change Initiative (1992-2015) oraz Copernicus Climate Change Service (2016-2020). Charakteryzują się one rozdzielczością przestrzenną równą 300x300 metrów. Ze względu na bardzo znikome lub zerowe zmiany pokrycia terenu, Antarktyda została wyłączona z analizy.
folders = c("africa", "asia", "europe", "north_america", "oceania", "south_america")
esa_files <- list(paste0("ESACCI-LC-L4-LCCS-Map-300m-P1Y-", 1992:2020, "-v2.0.7"))
for (i in 1:length(esa_files)) {
file_name <- paste0("esa_cci/", esa_files[i], ".tif")
rast <- rast(file_name)
year <- regmatches(esa_files[[i]], regexpr("\\d{4}", esa_files[[i]]))
assign(paste0("LC_", year), rast)
}1.2 Wczytanie danych wektorowych
Aby w pełni poprawnie obliczyć zmiany powierzchniowe, posłużono się danymi pochodzącymi z Equi7Grid - projektu odpowiadającego za utworzenie wiernoodległościowych systemów odniesienia dla wszystkich 7 kontynentów. Indywidualne siatki odwzorowań mają na celu zwiększenie efektywności przetwarzania wysokorozdzielczych danych przestrzennych.
for (i in 1:6) {
files = list.files(paste0(folders[i]), pattern = "\\GEOG_LAND.shp$", full.names = TRUE)
continents = vect(files)
continents_proj = project(continents, "+proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
assign(folders[i], continents)
}
cont_layers = list(africa, asia, europe, north_america, oceania, south_america)
layers_crs <- c("+proj=aeqd +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298 +datum=WGS84 +units=m +no_defs", # africa
"+proj=aeqd +lat_0=47 +lon_0=94 +x_0=4340913.84808 +y_0=4812712.92347 +datum=WGS84 +units=m +no_defs", # asia
"+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +datum=WGS84 +units=m +no_defs", # europe
"+proj=aeqd +lat_0=52 +lon_0=-97.5 +x_0=8264722.17686 +y_0=4867518.35323 +datum=WGS84 +units=m +no_defs", # north america
"+proj=aeqd +lat_0=-19.5 +lon_0=131.5 +x_0=6988408.5356 +y_0=7654884.53733 +datum=WGS84 +units=m +no_defs", # oceania
"+proj=aeqd +lat_0=-14 +lon_0=-60.5 +x_0=7257179.23559 +y_0=5592024.44605 +datum=WGS84 +units=m +no_defs") # south america
year <- 1992:2020
LC_list <- list(paste0('LC_', 1992:2020))2. Przycięcie oraz przekształcenie danych rastrowych
Założeniem przygotowania danych było uzyskanie warstw rastrowych pokrycia terenu każdego kontynentu dla każdego roku. Głównym etapem przygotowania danych do analizy było zatem przycięcie rastrów do obszarów poszczególnych kontynentów oraz zmiana układów współrzędnych z WGS84 na indywidualne odwzorowania dostarczone przez Equi7Grid.
Z powodu długotrwałego przekształcenia danych przy użyciu pakietów raster oraz terra, zdecydowano się wykorzystać w tym celu gdalUtils, który znacznie zmniejszył czas obliczeń. Etap ten był najdłużej trwającym etapem przygotowywania danych i zajął łącznie kilkadziesiąt godzin przetwarzania, głównie ze względu na problematyczny system odwzorowania dla Oceanii. W celu usprawnienia obliczeń zdecydowano się na podzielenie rastra Oceanii na dwie części, według równoleżnika 180°, co znacznie zmniejszyło czas transformacji układu z ponad 5 godzin do ok. 40 minut dla tego kontynentu.
W pracy magisterskiej zapewne warto umieścić tabelę przedstawiającą czas, jaki zajęły poszczególne procesy.
transf_oceania <- function(){
## ------------ first part -----------
rast_1 = crop(x, extent(70, 180, -60, 30))
raster_cropped_1 <- paste0("projected/", folders[j], "_temp1", year, ".tif")
writeRaster(rast_1, filename = raster_cropped_1, datatype = "INT1U")
gdalwarp(srcfile = raster_cropped_1,
dstfile = paste0("projected/", folders[j], "_proj0", year, ".tif"),
t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
gdal_translate(paste0("projected/", folders[j], "_proj0", year, ".tif"), paste0("projected/", folders[j], "_proj_1_", year, ".tif"),
of = "COG", co = "compress=lzw")
## ---------- second part ---------
rast_2 = crop(x, extent(-180, -125, -50, 10))
raster_cropped_2 <- paste0("projected/", folders[j], "_temp2", year, ".tif")
writeRaster(rast_2, filename = raster_cropped_2, datatype = "INT1U")
gdalwarp(srcfile = raster_cropped_2,
dstfile = paste0("projected/", folders[j], "_proj1", year, ".tif"),
t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
gdal_translate(paste0("projected/", folders[j], "_proj1", year, ".tif"), paste0("projected/", folders[j], "_proj_2_", year, ".tif"),
of = "COG", co = "compress=lzw")
}
transf_cont <- function(){
raster_cropped <- paste0("projected/", folders[j], "_temp.tif")
writeRaster(x, filename = raster_cropped, datatype = "INT1U")
gdalwarp(srcfile = raster_cropped,
dstfile = paste0("projected/", folders[j], "_proj0", year, ".tif"),
t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
gdal_translate(paste0("projected/", folders[j], "_proj0", year, ".tif"), paste0("projected/", folders[j], "_proj_", year, ".tif"),
of = "COG", co = "compress=lzw")
}Przygotowano dwie osobne funkcje transf_oceania oraz transf_cont, odpowiadające za zmianę układu współrzędnych przyciętych rastrów. Funkcje te wykorzystano w poniższej pętli dla wszystkich plików rastrowych oraz wszystkich lat.
for(i in LC_list) {
for (j in 1:length(cont_layers)) {
x = crop(i, cont_layers[j])
x = terra::mask(x, vect(cont_layers[j]))
if(j == 5){
transform_oceania()
file.remove(paste0("projected/", folders[j], "_proj0", year, ".tif"))
file.remove(paste0("projected/", folders[j], "_proj1", year, ".tif"))
file.remove(raster_cropped_1)
file.remove(raster_cropped_2)
oceania_1 <- rast(paste0("projected/", folders[j], "_proj_1_", year, ".tif"))
oceania_2 <- rast(paste0("projected/", folders[j], "_proj_2_", year, ".tif"))
merged <- vrt(c(oceania_1, oceania_2), paste0("projected/", folders[j], "_proj_", year, ".tif"))
file.remove(paste0("projected/", folders[j], "_proj_1_", year, ".tif"))
file.remove(paste0("projected/", folders[j], "_proj_2_", year, ".tif"))
} else{
transf_cont()
file.remove(paste0("projected/", folders[j], "_proj0", year, ".tif"))
file.remove(raster_cropped)
}
}
year <<- year + 1
}3. Reklasyfikacja
Obliczenie zmian złożoności przestrzennej zawartej w temacie pracy magisterskiej wymagało trzykrotnej reklasyfikacji danych na dwie klasy: obszary leśne oraz pozostałe, obszary zurbanizowane oraz pozostałe, a także obszary rolnicze oraz pozostałe. Zabieg ten pozwolił na szczegółowe zbadanie zmian struktury przestrzennej lokalnych krajobrazów oraz głównych trendów, m. in.: wylesiania, ekspansji rolnictwa oraz urbanizacji.
Reklasyfikację przeprowadzono w oparciu o zgrupowane klasy pokrycia terenu, proponowane przez IPCC (str. 30) w celu uniwersalnego porównania zmian 6 głównych klas pokrycia terenu, a także w oparciu o dokumentację klas ESA CCI.
m_forest <- c(0, 49, 0,
49, 101, 1,
101, 159, 0,
159, 171, 1,
171, 221, 0
)
matrix_forest <- matrix(m_forest, ncol = 3, byrow = TRUE)
years <- 1992:2020
folders = c("africa", "asia", "europe", "north_america", "oceania", "south_america")
for (i in 1:length(years)) {
files <- list.files(paste0(years[i]), pattern = ".tif$", full.names = TRUE)
for (j in 1:length(files)) {
continent <- rast(files[j])
raster_forest <- classify(continent, matrix_forest, include.lowest = TRUE, brackets=FALSE)
writeRaster(raster_forest, paste0("reclassified/", folders[j], "_forest_", years[i], ".tif"), datatype = "INT1U", overwrite = TRUE)
}
}Dane otrzymane po reklasyfikacji są zbiorem rastrów w podziale na lata, kontynenty oraz klasy pokrycia terenu.
library(terra)
afr_1992 <- rast("reclassified/1992/south_america_forest_1992.tif")
afr_2020 <- rast("reclassified/2020/south_america_forest_2020.tif")
plot(c(afr_1992, afr_2020), main = c("1992", "2020"), col = c("#444444", "#468b69"))
Przebieg analizy
Dalsze etapy pracy obejmują wykorzystanie zreklasyfikowanych danych do obliczenia metryk złożoności przestrzennej oraz udziału procentowego poszczególnych klas pokrycia terenu, w podziale na kontynenty oraz lata.
Obliczenie metryk
Wartości entropii brzegowej oraz względnej informacji wzajemnej otrzymano za pomocą kluczowych funkcji z pakietu motif. W celu ujednolicenia wyników, a także metod obliczenia statystyk i wykonania map, utworzono obszerną ramkę danych zawierającą niezbędne wyniki oraz informacje. Metryki wyliczono dla ponad 160 tys. lokalnych krajobrazów dla całego świata, każdy o rozdzielczości przestrzennej 30x30 km (100 razy mniejszej od rozdzielczości danych surowych). Otrzymana ramka metrics zawiera informacje o ID lokalnego krajobrazu, wartości entropii brzegowej, informacji wzajemnej oraz względnej informacji wzajemnej, a także o przyjętej rozdzielczości przestrzennej oraz procentowym udziale wartości NA, klasy pokrycia terenu i obszarów pozostałych. Dla każdego z nich dołączono również jego geometrię, uzyskując obszerny obiekt sf.
calc_relmutinf = function(mutinf, ent){
ifelse(mutinf == 0, 1, mutinf / ent)
}
add_it_metrics = function(x, year, res){
x %>%
mutate(ent = map_dbl(signature, comat:::rcpp_ent)) %>%
mutate(mutinf = map_dbl(signature, comat:::rcpp_mutinf)) %>%
mutate(relmutinf = calc_relmutinf(mutinf, ent)) %>%
mutate(res = res) %>%
select(-signature) %>%
set_names(c("id", "naprop", "ent", "mutinf", "relmutinf", "res"))
}
calculate_all_metrics = function(year, continent, class){
input_data = read_stars(file)
window = 100
res = 30
input_signature = lsp_signature(input_data, type = "coma", window = window)
input_metrics = add_it_metrics(input_signature, year, res = res)
input_metrics <- lsp_add_sf(input_metrics)
input_proc = lsp_signature(input_data, type = "composition", window = window)
input_proc = lsp_restructure(input_proc)
input_proc_sf = lsp_add_sf(input_proc)
metrics <- cbind(input_metrics,
"prop_0" = input_proc$X1,
"prop_1" = input_proc$X2,
"year" = year,
"continent" = continent,
"class" = class)
names(metrics)[names(metrics) == "naprop"] <- "prop_NA"
names(metrics)[names(metrics) == "prop_0"] <- "prop_0"
names(metrics)[names(metrics) == "prop_1"] <- "prop_1"
metrics
}W celu jak najdokładniejszej identyfikacji lokalnych krajobrazów na późniejszych etapach analizy, na podstawie nazewnictwa rastrów dodano informację o kontynencie, roku oraz reprezentowanej klasie pokrycia terenu. Poniższą pętlę powtórzono dla każdej klasy pokrycia terenu, co łącznie przełożyło się na ponad 35 godzin przetwarzania danych.
years <- 1992:2020
metrics <- data.frame()
for (i in 1:length(years)) {
files <- list.files(paste0(years[i]), pattern = "_forest_.*\\.tif$", full.names = TRUE)
for (file in files) {
file_name <- basename(file)
num_tokens <- length(unlist(strsplit(file_name, "_")))
if (num_tokens == 4) {
continent <- gsub("^([^_]+_[^_]+)_.*", "\\1", file_name)
class <- gsub("^.*?_[^_]+_([^_]+)_.*", "\\1", file_name)
}
else if (num_tokens == 3) {
continent <- gsub("^(.*?)_.*", "\\1", file_name)
class <- gsub("^.*?_(.*?)_.*", "\\1", file_name)
}
year <- as.integer(gsub("^.*?_.*_(\\d{4}).*", "\\1", file_name))
df <- calculate_all_metrics(year, continent, class)
df <- as.data.frame(df)
metrics <- rbind(metrics, df)
}
}metrics <- sample_n(metrics, 10)
subset(metrics, select = -geometry)| id | prop_NA | ent | mutinf | relmutinf | res | prop_0 | prop_1 | year | continent | class |
|---|---|---|---|---|---|---|---|---|---|---|
| 23038 | 0.009 | 0.0000 | 0.0000 | 1.0000 | 30 | 1.0000 | 0.0000 | 1992 | africa | forest |
| 80327 | 0.000 | 0.0000 | 0.0000 | 1.0000 | 30 | 1.0000 | 0.0000 | 2013 | africa | urban |
| 60486 | 0.000 | 0.3325 | 0.2172 | 0.6531 | 30 | 0.9393 | 0.0607 | 2001 | africa | agro |
| 53615 | 0.000 | 0.2559 | 0.0851 | 0.3325 | 30 | 0.0432 | 0.9568 | 1993 | africa | forest |
| 439842 | 0.000 | 0.6787 | 0.4010 | 0.5908 | 30 | 0.8204 | 0.1796 | 2011 | asia | forest |
| 331751 | 0.000 | 0.1012 | 0.0380 | 0.3751 | 30 | 0.9868 | 0.0132 | 2007 | asia | agro |
| 268588 | 0.000 | 0.0000 | 0.0000 | 1.0000 | 30 | 1.0000 | 0.0000 | 2008 | north_america | agro |
| 359020 | 0.000 | 0.0000 | 0.0000 | 1.0000 | 30 | 1.0000 | 0.0000 | 2000 | asia | urban |
| 369952 | 0.000 | 0.0000 | 0.0000 | 1.0000 | 30 | 1.0000 | 0.0000 | 2018 | asia | urban |
| 125687 | 0.000 | 0.0291 | 0.0162 | 0.5586 | 30 | 0.9969 | 0.0031 | 2018 | oceania | urban |
Grupowanie metryk
W celu obliczenia zmian wartości metryk, z otrzymanego obiektu sf metrics wydzielono lata 1992 i 2020 oraz pogrupowano według kontynentu, klasy i geometrii, aby dodać kolumnę z unikalnym ID każdego lokalnego krajobrazu. Pogrupowane rekordy dla różnych lat przyjęły identyczne ID - to pozwoliło na obliczenie zmiany wartości entropii brzegowej oraz względnej informacji wzajemnej na podstawie par krajobrazów lokalnych. W celu ich późniejszej identyfikacji dodano informację o kontynencie i klasie pokrycia terenu, a także pierwotnym ID krajobrazu.
sf_result <- metrics %>%
group_by(continent, class, geometry) %>%
mutate(group_id = cur_group_id())
df_grouped <- sf_result %>%
group_by(group_id) %>%
summarize(ent_diff = diff(ent),
mutinf_diff = diff(mutinf),
relmutinf_diff = diff(relmutinf))
df_grouped$continent <- sf_result$continent[match(df_grouped$group_id, sf_result$group_id)]
df_grouped$class <- sf_result$class[match(df_grouped$group_id, sf_result$group_id)]
df_grouped$id <- sf_result$id[match(df_grouped$group_id, sf_result$group_id)]df_grouped <- st_drop_geometry(df_grouped)
df_grouped <- df_grouped[order(df_grouped$id),]
df_grouped <- df_grouped[1:9,]| group_id | ent_diff | mutinf_diff | relmutinf_diff | continent | class | id |
|---|---|---|---|---|---|---|
| 419470 | -0.0473 | -0.0103 | -0.0147 | south_america | agro | 175 |
| 440163 | 0.0691 | 0.0545 | 0.0465 | south_america | forest | 175 |
| 460856 | 0.1993 | 0.0433 | -0.0690 | south_america | urban | 175 |
| 419471 | 0.0126 | 0.0008 | -0.0043 | south_america | agro | 788 |
| 440164 | -0.0002 | -0.0074 | -0.0073 | south_america | forest | 788 |
| 460857 | 0.0162 | 0.0090 | 0.4089 | south_america | urban | 788 |
| 419472 | 0.0425 | 0.0619 | 0.0625 | south_america | agro | 1088 |
| 440165 | 0.0224 | 0.0693 | 0.0641 | south_america | forest | 1088 |
| 460858 | 0.0852 | 0.0380 | 0.0846 | south_america | urban | 1088 |
Obliczenie powierzchni klas
W celu otrzymania szczegółowych zmian powierzchni, zreklasyfikowane rastry zamieniono na ramki danych, dodając do każdego piksela informacje identyfikacyjne. Ramkę tę następnie pogrupowano według kontynentu i roku, aby zsumować liczbę występowań pikseli i obliczyć udział procentowy każdej klasy.
years <- 1992:2020
result_df <- data.frame()
for (i in 1:length(years)) {
files <- list.files(paste0(years[i]), pattern = ".tif$", full.names = TRUE)
for (file in files) {
raster_data <- rast(file)
freq_data <- as.data.frame(freq(raster_data))
file_name <- basename(file)
num_tokens <- length(unlist(strsplit(file_name, "_")))
if (num_tokens == 4) {
freq_data$continent <- gsub("^([^_]+_[^_]+)_.*", "\\1", file_name)
freq_data$class <- ifelse(freq_data$value == 0, "other", gsub("^.*?_[^_]+_([^_]+)_.*", "\\1", file_name))
}
else if (num_tokens == 3) {
freq_data$continent <- gsub("^(.*?)_.*", "\\1", file_name)
freq_data$class <- ifelse(freq_data$value == 0, "other", gsub("^.*?_(.*?)_.*", "\\1", file_name))
}
freq_data$year <- as.integer(gsub("^.*?_.*_(\\d{4}).*", "\\1", file_name))
total_cells <- sum(freq_data$count)
freq_data <- freq_data %>%
group_by(continent, year) %>%
mutate(percentage = round((count / total_cells) * 100, 2)) %>%
ungroup()
result_df <- bind_rows(result_df, freq_data)
}
}Udział procentowy obszarów pozostałych był nieprawdziwy, ze względu na potrójne wystąpienia klasy “pozostałe” dla każdego roku i każdego kontynentu. W celu poprawnego obliczenia udziału procentowego, ramkę danych pogrupowano oraz zsumowano udział procentowy do 100, względem pozostałych prawdziwych klas.
result_df_updated <- result_df %>%
group_by(year, continent, class) %>%
mutate(total_percentage = sum(percentage)) %>%
group_by(year, continent) %>%
mutate(percentage = ifelse(class == "other", 100 - sum(percentage[class != "other"]), percentage)) %>%
ungroup() %>%
select(-total_percentage) %>%
distinct(year, continent, class, .keep_all = TRUE)Rzeczywista liczba pikseli klasy “pozostałe” została obliczona na podstawie różnicy zsumowanej liczby pikseli wszystkich klas oraz trzech klas głównych (obszary leśne, zurbanizowane, rolnicze). W celu prezentacji statystyk dla całego świata, dane te pogrupowano według roku i kontynentu, dodając informacje o zsumowanej liczbie pikseli i udziale procentowym.
continents <- c("africa", "asia", "europe", "north_america", "oceania", "south_america")
files <- list.files(pattern = "_.*_\\d{4}.tif$", full.names = TRUE)
files <- files[seq(1, length(files), by = 3)]
cont_proc <- c()
for (file in files) {
raster_data <- rast(file)
freq_data <- as.data.frame(freq(raster_data))
total_cells <- sum(freq_data$count)
cont_proc <- c(cont_proc, total_cells)
}
result_df_updated <- result_df_updated %>%
group_by(continent, year) %>%
mutate(count = ifelse(class == "other", cont_proc[match(continent, c("africa", "asia", "europe", "north_america", "oceania", "south_america"))] - sum(count[class %in% c("agro", "forest", "urban")]), count))
world_aggregated <- result_df_updated %>%
filter(continent != "world") %>%
group_by(year, class) %>%
summarise(count = sum(count),
percentage = round((count/sum(cont_proc)) * 100, 2),
continent = "world")
result_df_updated <- bind_rows(result_df_updated, world_aggregated)Taki zestaw danych pozwolił na obliczenie zmian powierzchni poszczególnych klas dla każdego kontynentu dla całego badanego okresu.
result_df_updated <- result_df_updated[1:10,]
result_df_updated| count | continent | year | class | percentage |
|---|---|---|---|---|
| 241287787 | africa | 1992 | other | 62.84 |
| 56765277 | africa | 1992 | agro | 14.78 |
| 85626299 | africa | 1992 | forest | 22.30 |
| 312359 | africa | 1992 | urban | 0.08 |
| 198747887 | asia | 1992 | other | 45.07 |
| 97441833 | asia | 1992 | agro | 22.10 |
| 143880631 | asia | 1992 | forest | 32.63 |
| 880366 | asia | 1992 | urban | 0.20 |
| 23817815 | europe | 1992 | other | 21.18 |
| 47150640 | europe | 1992 | agro | 41.93 |
Wyniki
Globalne zmiany powierzchniowe
Obliczenie zmian udziału powierzchni może ukazać nie tylko trendy poszczególnych klas pokrycia terenu dla każdego kontynentu, ale także najważniejsze momenty w historii, które dany trend zaburzały lub odwracały.
Długotrwały, stabilny trend wzrostowy obszarów leśnych jest zauważalny jedynie w Afryce. Od kilku lat występuje również w Ameryce Południowej. Oceania jako jedyny kontynent odznacza się trendem spadkowym od roku 1992.
Największa zmiana powierzchni obszarów leśnych, sięgająca ponad 500 tys. km2, nastąpiła w Ameryce Południowej. Zmiana ta jest związana przede wszystkim z długotrwałą wycinką Puszczy Amazońskiej. Największym jednorocznym spadkiem powierzchni charakteryzuje się rok 2004, w którym nastąpiło zmniejszenie obszarów leśnych w Ameryce Południowej o ponad 120 tys. km2. Zmiany powierzchni obszarów leśnych dla pozostałych kontynentów są znacznie mniejsze i w ciągu 28 badanych lat zmieniały swoje trendy.
Globalne zmiany złożoności przestrzennej
Obliczone metryki posłużyły do ukazania stanu złożoności przestrzennej, jaka występowała na całym świecie w każdym z badanych lat. Wizualnie zmiany te są zauważalne głównie w dużych odstępach czasowych, dlatego poniżej przedstawiono wartości metryk jedynie w dwóch skrajnych latach: 1992 i 2020 roku.
Dane w obiekcie metrics nie posiadały przypisanego układu współrzędnych, zatem w celu ich poprawnego wyświetlenia należało przypisać każdemu kontynentowi jego indywidualne odwzorowanie. Następnym etapem była rasteryzacja danych ze względu na niepoprawne wyświetlanie obiektów w formie poligonowej oraz zmiana układu współrzędnych na wspólny, docelowy układ (tutaj: WGS84).
metrics_map <- function(class, year, metric){
layers_crs <- c("+proj=aeqd +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298 +datum=WGS84 +units=m +no_defs", # africa
"+proj=aeqd +lat_0=47 +lon_0=94 +x_0=4340913.84808 +y_0=4812712.92347 +datum=WGS84 +units=m +no_defs", # asia
"+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +datum=WGS84 +units=m +no_defs", # europe
"+proj=aeqd +lat_0=52 +lon_0=-97.5 +x_0=8264722.17686 +y_0=4867518.35323 +datum=WGS84 +units=m +no_defs", # north america
"+proj=aeqd +lat_0=-19.5 +lon_0=131.5 +x_0=6988408.5356 +y_0=7654884.53733 +datum=WGS84 +units=m +no_defs", # oceania
"+proj=aeqd +lat_0=-14 +lon_0=-60.5 +x_0=7257179.23559 +y_0=5592024.44605 +datum=WGS84 +units=m +no_defs") # south america
cont_layers = list("africa", "asia", "europe", "north_america", "oceania", "south_america")
df <- c()
for (i in 1:6){
data <- df_sf[df_sf$continent == cont_layers[i],]
data <- st_set_crs(data, layers_crs[i])
df[[i]] <- data
}
rast_list <- list()
for (j in 1:length(df)){
rast <- df[[j]]
empty_raster <- raster(rast , res = 30000)
rast <- terra::rasterize(rast, empty_raster, field = metric)
rast <- projectRaster(rast , crs = "+proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs", res = 30000)
origin(rast) <- 0
rast_list[[j]] <- rast
}
merged <- do.call(merge, rast_list)
assign(paste0(metric, '_', class, '_', year), merged, envir = .GlobalEnv)
}
metrics_map("forest", "2020", "ent")
projected <- projectRaster(ent_forest_1992, crs = 4326)Entropia brzegowa
Rozkład przestrzenny wartości entropii brzegowej dla obu lat jest wizualnie bardzo podobny. Największa różnorodność krajobrazu leśnego występuje na większości obszarów Europy, a także w południowo-wschodniej Azji oraz na terenach tajgi w północnej Kanadzie. Rozkład wartości uległ zmianie - nastąpił lekki przyrost wszystkich obszarów z wartościami entropii powyżej 0, co oznacza zwiększenie zróżnicowania struktury krajobrazu.
Najbardziej widoczną zmianą średniej wartości entropii w ciągu 28 lat charakteryzuje się Ameryka Południowa., która zanotowała znaczący wzrost wartości metryki. Pozostałe kontynenty cechuje brak istotnej zmienności zróżnicowania przestrzennego krajobrazów. W ciągu ostatnich kilku lat wystąpił zauważalny trend wzrostowy entropii brzegowej dla każdego z kontynentów.
Największym wzrostem entropii brzegowej między 1992 a 2020 rokiem odznacza się Amazonia. Znaczące spadki są zauważalne również na obszernych terenach leśnych we wschodniej Rosji, a także północnej Nikaragui i południowo-zachodniej Kanadzie. Obszary charakteryzujące się wzrostem entropii brzegowej są dwukrotnie liczniejsze od obszarów, dla których wartość metryki zmalała. Największe spadki zróżnicowania przestrzennego krajobrazu odnotowano na wielu obszarach Ameryki Południowej oraz Afryki Środkowej, a także w Birmie oraz stanie Minnesota w USA.
Wszystkie kontynenty odznaczają się znaczną przewagą obszarów ze wzrostem entropii brzegowej od obszarów, gdzie uległa ona spadkowi. Największa liczba lokalnych krajobrazów, dla których entropia brzegowa uległa znacznemu wzrostowi, występuje w Ameryce Południowej. Afryka odznacza się największym udziałem procentowym obszarów, które nie uległy zmianie, Europa natomiast - najmniejszym.
Względna informacja wzajemna
Przestrzenne zróżnicowanie wartości względnej informacji wzajemnej pomiędzy rokiem 1992 i 2020 jest dość widoczne. Dla obszaru całego świata nastąpił spadek wartości metryki, a więc zmniejszyła się jednolitość struktury przestrzennej lokalnych krajobrazów. Największy spadek wartości względnej informacji wzajemnej jest najbardziej widoczny na obszarze Amazonii, największy wzrost natomiast - w północnych Indiach.
Najbardziej widoczna zmianą średniej wartości względnej informacji wzajemnej wystąpiła w Ameryce Południowej w tym samym roku (2004), w którym jej średnie wartości entropii brzegowej znacznie wzrosły. Kontynentem znacznie odstającym od reszty jest Afryka, ze średnią wartością względnej informacji wzajemnej powyżej 0,6. W ciągu ostatnich kilku lat, w kontraście do rosnącej entropii, względna informacja wzajemna ulega zmniejszeniu na każdym z kontynentów.
Największym spadkiem względnej informacji wzajemnej między 1992 a 2020 rokiem odznacza się Amazonia. Znaczące spadki wystąpiły również w Argentynie, środkowej części USA oraz na granicy Sudanu Południowego i Demokratycznej Republiki Konga. Łącznie na świecie więcej obszarów odznaczyło się spadkiem względnej informacji wzajemnej, aniżeli wzrostem. Największy wzrost obszarów, które w ciągu 28 lat stały się bardziej jednolite, nastąpił w północno-wschodnich Indiach, wschodnich Chinach oraz w południowo-zachodniej części Ameryki Południowej.
Europa jest jedynym kontynentem, którego większość obszarów odznacza się wzrostem względnej informacji wzajemnej. Pomimo dużego spadku jednolitości struktury przestrzennej na większości obszarów w Ameryce Południowej, wzrost względnej informacji wzajemnej również wystąpił na znacznej liczbie obszarów i jest relatywnie niewiele mniejszy.
Lokalne zmiany złożoności przestrzennej
Badanie zmian złożoności przestrzennej na poziomie lokalnym może znacznie przybliżyć przyczyny występujących zdarzeń i trendów, które wpłynęły na zwiększenie zróżnicowania struktury przestrzennej krajobrazów lub zwiększenia jej jednolitości.
Entropia brzegowa
Najwyższy spadek entropii brzegowej
Obszary z najwyższym spadkiem entropii brzegowej występują w południowej oraz wschodniej części Puszczy Amazońskiej w Brazylii. W ciągu 28 lat ponad 6 mln km2 obszarów przedstawionych na poniższej rycinie uległo zwiększonemu rozdrobnieniu krajobrazu. Tylko na ponad 600 tys km2 obszarach entropia brzegowa uległa lekkiemu spadkowi.
Największy wzrost entropii brzegowej
Znaczący wzrost jednolitości struktury krajobrazu wystąpił na granicy Sudanu Południowego oraz Demokratycznej Republiki Konga w Afryce. Zmiany te dotyczą prawie 300 tys. km2 obszarów, dla których nastąpił spadek entropii brzegowej - wystąpiło zalesienie wielu obszarów, które przedstawiono poniżej. Obszar ten wyróżnia się na tle terenów sąsiadujących, na których entropia brzegowa uległa wzrostowi.
Względna informacja wzajemna
Najwyższy spadek względnej informacji wzajemnej
Oprócz najwyższego spadku zróżnicowania struktury krajobrazu, najwyższy spadek jednolitości struktury również występuje w Puszczy Amazońskiej. Obszary te są natomiast rozproszone po całej Amazonii, bez znacznego skupiska w Brazylii. W ciągu 28 lat ponad 5,5 mln km2 zanotowało spadek względnej informacji wzajemnej, co związane jest z długoletnią wycinką postępującą na obszarze Amazonii. W sąsiedztwie występuje również ok. 2 mln km2 obszarów, na których nastąpił lekki wzrost metryki.
Największy wzrost względnej informacji wzajemnej
Obszary ze znaczącym wzrostem jednolitości struktury krajobrazu znajdują się na północy Indii i są związane z ubytkiem obszarów leśnych. Główne zmiany dotyczą ok. 150 tys. km2, na których zwiększyła się wartość względnej informacji wzajemnej.
Identyfikacja obszarów o podobnej złożoności przestrzennej
…
Podsumowanie
W ciągu 28 badanych lat na zdecydowanej większości obszarów leśnych na świecie nastąpiła zmiana złożoności przestrzennej. Łącznie na świecie znajduje się znacznie więcej obszarów, na których różnorodność struktury krajobrazu uległa zwiększeniu, niż obszarów, gdzie wystąpił ich wzrost jednolitości. Świadczy to o postępującym rozdrabnianiu krajobrazu leśnego, co w ostatnich latach postępuje najmocniej. Głównym obszarem największych zmian w różnorodności oraz jednolitości struktury krajobrazu jest Puszcza Amazońska, gdzie od wielu lat występuje wycinka drzewostanu leśnego, która w ostatnich latach spowolniła i jest rekompensowana przez zalesianie.